Nonequilibrium phase transition in the sedimentation of reproducing particles 
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We study numerically and analytically the dynamics of a sedimenting suspension of active, re- 
producing particles, such as growing bacteria in a gravitational field. In steady state we find a 
non-equilibrium phase transition between a 'sedimentation' regime, analogous to the sedimentation 
equilibrium of passive colloids, and a 'uniform' regime, in which the particle density is constant in all 
but the top and bottom of the sample. We discuss the importance of fluctuations in particle density 
in locating the phase transition point, and report the kinetics of sedimentation at early times. 
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About a century ago, Einstein showed theoretically 
and Pcrrin demonstrated experimentally that in a dilute 
colloidal suspension, the particle density, n, as a func- 
tion of height, z, is given by the barometric distribution: 
n(z) = n(0) exp(— z/zq), where zq is the sedimentation 
height [H. This distribution results from a subtle in- 
terplay between thermal diffusion, hydrodynamics and 
gravity. Diffusion and hydrodynamics are related via 
the Stokes-Einstein formula, a form of the fluctuation- 
dissipation theorem for equilibrium systems, while bal- 
ancing gravity and Brownian motion gives zq = D/vs, 
where D is the particles' diffusion coefficient and Vs their 
sedimentation speed. The barometric distribution ap- 
plies when the suspension is so dilute that interparticle 
potential interactions (excluded volume. Coulomb, etc.) 
can be neglected. Modern colloid physics has focussed 
on the behavior of concentrated suspensions j3l . 

An equally interesting avenue to explore is that of ac- 
tive particles ( APs) |3l . Specifically, we consider APs able 
to propel themselves in such a way that their long-time 
motion is diffusive, i.e. each particle's mean-squared dis- 
placement from an initial position, (r^(i)), increases lin- 
early with time, t, so that (in three dimensions) (r^) = 
6-Doffi, where Deff is an effective diffusion coefScient. 
The swim-tumble-swim motion of an Escherichia coli 
bacterium (~ 2yumxl/im splieroclinder, average density 
Pb = 1.08 g/cvo?) is an example l4|], for which experi- 
ments give Doff ~ C(10^/xm^s~^) [5|. An equivalent pas- 
sive colloid has D ~ 0.5 /im^s~^ at 300 K, so that an E. 
coli functions at an effective temperature of 0(10^ K): it 
is far from equilibrium. Mimicking bacteria, we assume 
our APs can also 'reproduce' and 'die': their number 
density may therefore change with time. Both motility 
and reproduction require energy intake, although this is 
nowhere explicit in what follows. 

In this Letter, we study the behavior of (effectively) 
diffusing and reproducing, non-interacting APs in a grav- 
itational field. This may model, for instance, a dilute 
suspension of motile E. coli that are growing but not re- 
sponding to chemical gradients (i.e. non-chemotactic) . 
Such a system is, arguably, the simplest example of 'ac- 
tive soft matter'. Does this paradigmatic AP system dif- 



fer significantly from its passive counterpart (dilute col- 
loidal sedimentation equilibrium)? And if so, how? 

We first report stochastic simulations of the dynamics 
of diffusing and reproducing APs in a gravitational field. 
We then interpret these results by analysing a continuum 
equation describing the evolution of a density profile of 
a dilute sedimenting AP fluid. We find the steady state 
profile, as well as the dynamic pathway leading to it. 
Even a mean-field description reveals a much richer phe- 
nomenology for APs than passive colloids. We find a 
non-equilibrium phase transition between a 'sedimenta- 
tion' regime with exponential profile, and another regime 
showing essentially constant density in the bulk of the 
suspension. Using realistic parameter values, we predict 
that one may switch between the two phases by modify- 
ing the growth rate of a real system of bacteria. Close to 
this transition, there exist novel 'sedimentation bands' 
in which a region of uniform AP density coexists with 
an AP-depleted region. This may usefully be compared 
to the phenomenon of shear banding. We also analyse 
the role of noise, and show that its presence shifts the 
transition point. Finally, we discuss generalisations of 
our equations, and how our results may relate to real 
bacterial suspensions. 

We use a stochastic algorithm to simulate the coupled 
biased diffusion and reproduction/death of APs. We con- 
sider a column of sedimenting APs as a discrete lattice of 
sites z = 1 , . . . L with the number of AP occupying each 
lattice site specified as Uiit). Gravity acts towards i = 0. 
At each time step t t-\-/S.t the array of occupation num- 
bers is updated according to a 'multiply' step or a 'move' 
step chosen with probabilities w/{\-\-w) or 1/(1 -I- w) re- 



spectively where w 



^ [l + Ef=i"f/(PoE^"^W)J is 
the ratio of the total rate of reproduction/death per par- 
ticle to the total rate of moving per particle, and a, po 
are parameters to be discussed below 0] . 

In a 'move' update each AP moves independently up or 
down with probability p or 1 — p. In a 'multiply' update, 
at each site i each particle is replaced by two particles 
with probability po/{po + n'i{t)) or removed with proba- 
bility ni{t) / {po + ni{t)) . We impose no flux, or reflecting, 
boundary conditions at the top and the bottom of the 
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container (i = and i = L). 

If Az, At represent respectively the spatial and tempo- 
ral steps, then the continuum limit of our Markov process 
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leads to a diffusion constant D 
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and sedimentation 



velocity v = {1 — 2p) The parameter p controls the 
strength of the gravitational force; a controls the rate 
of reproduction/death and po gives the value of the oc- 
cupation in which reproduction and death are balanced. 
In a real bacterial suspension, a will be the growth rate 
(medium dependent, but ^ hour"^ for E. coli) and Pn is 



the saturation cell density 10*^ cells/cm"^ for E. coli[^). 

Note that the state where the lattice is devoid of APs 
is an absorbing state of the dynamics in the algorithm, 
which we call model I. We also considered a variation of 
the dynamics, which we refer to as model II, in which 
a depopulated site is allowed to be repopulated sponta- 
neously (during the reproduction/death move). 

Our simulations suggest that there is a phase transi- 
tion between two different regimes. If the APs reproduce 
slowly (a small), we obtain a steady-state density profile 
which decays sharply with z (Fig. la). We call this the 
'sedimentation regime' in analogy with the similar be- 
haviour of passive colloids under gravity. If the growth 
rate exceeds a threshold, ac, the steady state is one with a 
uniform density throughout the bulk of the sample, with 
depleted and enriched layers in the top and the bottom 
of the container (Fig. lb). In model II, we find that the 
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critical value is close to ac whereas in model I, the 

presence of an absorbing state widens the sedimentation 
regime, shifting ac to slightly larger values. 

To understand the simulations, we take a continuum 
approach and coarse grain our discrete model by defin- 
ing a 1-dimcnsional profile of AP density, p{z, t) (gravity 
acts along the negative z axis). To make contact with 
the numerical simulations, and on general grounds, wc 
may assume that the density profile obeys the following 
dynamic equation of motion: 
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The first term represents diffusion, the second term repre- 
sents the sedimentation due to gravity and the third term 
represents reproduction/death. Here, D and v have the 
same meaning as in the numerical simulations. An im- 
portant dimensionless control parameter is 9 = vjyj Da. 
The last term in ([T]) represents noise, and ?7(t) is a white 
noise with unit variance, while r(p) is a function to 
be specified. We first consider the deterministic case, 
r(p) = which reduces to a mean-field description. 

Considering the steady state (t oo), if a = 0, 
we get back passive colloid sedimentation equilibrium 
p[z) oc exp(— wz/Z?) . However, for a > one obtains 
a nonlinear equation for which an exact solution is not 
available. Therefore we perform a perturbation expan- 
sion for small a using the first 30 terms of a series ex- 
pansion p(z) = X^n '^"/'^"H-^^)- The series converges for 
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FIG. 1: The average density, p, versus height, z, for a sus- 
pension of motile and reproducing active particles in a gravi- 
tational field. Simulation results (model I), with (a) d — 3.6, 
and (b) 6 — 1.4. In (c) we show diflferent numerical solutions 
of the noiseless continuum equation Eq. [T] while (d) shows 
the dependency of the sedimentation length, on the velocity 
(in units of V Da; the dashed line gives the colloidal value at 
a — 0, which is denoted by zq in the text). 



a < ac = jj^ and diverges otherwise. This corroborates 
the simulation results and suggests that indeed a phase 
transition occurs as the 9 ~ vj^J Da goes through 2. 

Eq. [T] with w = and no noise is the well-known 
Fisher-KPP equation and admits advancing waves with 
velocity Vyj = 2\fDa as solutions we can under- 
stand this phase transition as a competition between 
gravity forcing the bacteria downwards with velocity v 
and a travelling wave of proliferation which advances up- 
ward. For V > the sedimentation wave wins over the 
Fisher wave and leads to an exponential profile, while 
for V < the Fisher wave leads to a uniform bacterial 
density throughout the sample. However, Eq. [factually 
yields a transient travelling wave for v > only. The 
transition is reminiscent of transitions in interface depin- 
ning l9|; in a branching random walk with an absorbing 
wall llOl l , and of extinction transitions in inhomogeneous 
We note that a linearised version 
= and different boundary con- 
ditions was considered in Ref. [13| to describe microor- 



biological systems [U 
of our Eq. ^ with F 



ganisms advccted in a river and resultant extinction. We 
believe that our analysis should apply to this problem as 
well. We also stress that the transition would be washed 
away by translational invariance or periodic boundary 
conditions, which are usually considered (s. 11|. 

The existence of a non-equilibrium phase transition at 
9 = 9c = 2 is confirmed by numerical solution of the 
noiseless version of Eq. [1] using a standard finite dif- 
ference scheme. An example of a series of steady state 
solutions for different values of 9 is shown for a sample 
of size L = 20-\/D/a, in Fig. Ic. Increasing the sys- 
tem size, the segment of the sample at p = 1 in the 
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FIG. 2: Steady state diagram in the {v, L) plane for the noise- 
less version of Eq. [T] A density profile is classified as 'banded' 
if the cline stays in the bottom 75% of the sample. Typical 
concentration profiles in the various regimes are also shown, 
together with a color scale for the density. 
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FIG. 3: (a) Dynamics of the density profile in a system with 
L = 30 (in units of ^/Dja), and with 9 = 2.1. The times 
corresponding to each profile are (from left to right) 0.99, 2.8, 
17.5, 62.2, 121.4, 250 (in units of a'^ , the bottom 5% of the 
sample is cut for colour readability), (b) Plots of the time 
needed to get to steady state (in units of a~^), as a function 
of V (in units of \/Da), with a system with L = 20 (solid 
line), L = 28 (dashed line) and L — 40 (dot-dashed line). (L 
is in units of y^D/a.) 



uniform regime increases (data not shown) , analogous to 
equilibrium phase transitions. It is interesting to con- 
sider the behavior of the decay length of the exponential 
density profile in the sedimentation regime (an effective 
sedimentation length) as a function of \9 — 6c\. This sed- 
imentation length is akin to a scaling length in an equi- 
librium phase transition. We find that after correcting 
for a small L-dependent shift in 9c, the sedimentation 
length does not diverge at the transition, and is only at 
most ~ 20% larger than the corresponding sedimentation 
length with no growth (a = 0). If we focus on the steady 
state concentration value e.g. in the middle of the sam- 
ple, it switches abruptly, for L —f oo, from 0, for 9 > 9c, 
to 1, ior 9 < 9c. These obervations are consistent with a 
discontinous phase transition. 

For large but finite systems, we find evidence of an 
intriguing spontaneous banding of the sedimenting APs, 
which occurs very close to 9 = 9c. An example is shown 
in Fig. 2, in which two steady-state regions coexist in the 
bulk of the sample, one in which p is practically and 
another one in which it is ^ 1. This kink-like solution, 
which we call a 'sedimentation band', is similar to a cline, 
found in the population biology literature when solving 
a reaction-diffusion equation, similar to Eq. [Tl but with 
V = and a reaction term which is cubic in p [l3| . In our 
case, however, the existence of sedimentation bands is 
more surprising, as without advection the state at p = 
is unstable, and bands arise due to the vicinity of a phase 
transition. In this respect, our sedimentation banding is 
more akin to shear banding, which is obtained when some 
complex fluids such as liquid crystals and worm-like mi- 
celles in the isotropic or disordered phase, are subjected 
to a shear, slightly smaller than that needed to order 
them completely |14{ . Sedimentation bands appear in 
an increasingly small window of as L increases, and 
disappear in the thermodynamic (L oo) limit. 

Would it be possible to observe the transition we pre- 



dict in a real bacterial suspension? For E. coli in water, 
V ^ 0.1//ms~^ and is fixed, while D > 10^/xms~^ 
In rich, well aerated media maintained at the opti- 
mal temperature of 37"^ C, the population doubles every 
20 minutes, giving a ~ 10"'^ s~^, and 9 w 0.3. It is 
possible to culture the bacterium in what is known as 
'motility buffer', in which D is maintained, but growth 
essentially stops {a ^ 0), allowing the tuning of 9 from 
0.3 through 2 to arbitrarily large values, thus permitting 
the observation of our transition in principle. 

In this context, it is important to note that the path- 
way to steady state may be quite slow. This is partic- 
ularly true close to 9c{L), and in the region where sedi- 
mentation bands form. In the banding regime for large 
sample size L, wc also find that the behaviour of the part 
of the sample close to the top, or just after the boundary 
of the band, display non-monotonic behavior. The den- 
sity first increases, as if the systems transiently entered 
the uniform regime, to decay later on to reach equilibrium 
(see Fig. 3a) . The time scale needed to reach equilibrium, 
toq, is plotted in Fig. 3b as a function of the distance from 
the transition point. Larger systems take longer to equi- 
librate, while close to criticality we find that toq increases 
as a power law oi \9 — 9r\ with a ~ 1 above the tran- 
sition, consistent with [lO| and confirming the presence 
of a phase-transition at 0c = 2. These results applied to 
APs with a = 10~^ s~^ predict that close to the transi- 
tion, it may take up to several months for a column of 
10 cm height to reach steady state. 

Next, we discuss the role of noise, i.e. the case T{p) ^ 
in Eq. [TJ If p — is to be an absorbing state, we need to 
go beyond the case T{p) = Fq > 0, and consider instead 
a density-dependent amplitude F. The natural choice is 
r = Foy^, which is justified by a central limit argument 
relating the variance of the noise to the number of active 
particles[ll,[i3|. Increasing Fq favours large fluctuations 
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FIG. 4: The dependence of the critical velocity, Vc (in units 
of -/Da) on To (in units of {Daf^), for L = 200 ^/D/^. 



and may locally bias the system towards the absorbing 
state. We observe that as Tq is increased the sedimen- 
tation regime is enhanced at the expense of the growth 
one (see Fig. 4). This is consistent with a negative shift 
in the Fisher wave velocity due to noise [l^. The non- 
equilibrium phase boundary shown in Fig. 4 is found by 
locating the maxima of the order parameter fluctuations. 

Eq. [T]with u = describes directed percolation (DP) 
[fit , the generic university class of nonequilibrium tran- 
sitions from an absorbing state (p = 0) to a fluctuating 
state {p > 0). The dimensionless control parameter is in 
this case To/iDaf^'^. Thus Fig. 4 becomes an extension 
of the DP phase diagram to include a v axis. Intriguingly 
the transition we have studied for Tq is discontinuous 
whereas the DP transition at w = is continuous. There- 
fore one may speculate that there is a singular or even 
tricritical point along the critical curve Vc{^o)- This sce- 
nario might be similar to what occurs in the XY model, 
where switching the XY spin into a velocity can result in 
a discontinuous transition [20| . 

We have already pointed out that with typical values 
of V and D in a suspension of E. coli, it may be possible 
to 'tune' the growth constant a to bring the system from 
the sedimentation to the uniform regime according to our 
predictions. In reality, these experiments need stringent 
controls, e.g. to make sure that the bacteria are not 
engaged in any form of chemotaxis, which would ren- 
der D dependent on concentrations of chemical species 
(nutrient, oxygen, waste products, . . . ). The extra level 
of complexity introduced can be modelled by adding a 
chemotactic term to Eq. (TJ and coupling it to a reaction- 
diffusion equation, e.g. as in the Keller-Segel model [l9| . 
Moreover, we have shown that the time scales for reach- 
ing steady state can be long, reaching 0(10'^) in units 
of the inverse growth rate, a^^, in the vicinity of the 
transition. But in a bacterial culture, a itself is only 
approximately constant during what is known as the 'ex- 
ponential' growth phase, after which saturation in popu- 
lation density and then death follow. Thus, steady-state 
experiments at 6 ^ 6c are likely impractical. 

Clearly, our model also neglects hydrodynamic inter- 
actions 2]J, which may (for example) cause swimmers 



to attract. Thus, hydrodynamics could have highly non- 
trivial effects, e.g. concerning the approach to steady 
states. However, we believe that the qualitative features 
of the transition we have identified may survive, because 
these are due to the competition between gravity and 
Fisher wave fronts, which should be generic. Also, there 
are non-swimming bacteria, such as e.g. Sinorhizobium 
Melitoti, where hydrodynamics is less important. 

To summarize, we have shown that the physics of 
motile and 'reproducing' active particles in a gravita- 
tional field yields a surprisingly non-trivial phenomenol- 
ogy. We have found, by simulations and analytic work, 
that increasing the growth rate from zero, the system 
makes an abrupt transition from a sedimentation regime 
in which the density decays exponentially with the dis- 
tance from the bottom of the container, to a uniform 
growth-dominated regime in which the density is prac- 
tically constant spatially except for two boundary layers 
at the top and at the bottom to satisfy boundary condi- 
tions. The essential physics is that of a balancing between 
a downward gravitational flux, and an upward diffusion- 
growth flux. Using values appropriate to E. coli, we pre- 
dict that it may be possible to observe this transition in 
a real bacterial suspension. 
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